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Abstract 

Ecolab and Webworld are both models of evolution produced 
by adding evolution to ecological equations. They differ pri- 
marily in the form of the ecological equations. Both models 
are self-organised to a state where extinctions balance specia- 
tions. However, Ecolab shows evidence of this self-organised 
state being critical, whereas Webworld does not. This paper 
examines the self-organised states of these two models and 
suggest the likely cause of the difference. Also the lifetime 
distribution for a mean field version of Ecolab is computed, 
showing that the fat tail of the distribution is due to coevolu- 
tionary adaption of the species. 

Introduction 

In models of evolving ecologies, a "drip feed" of mutated 
species are added to a simulation of ecological dynamics. 
As new species are incorporated into the ecology, they cre- 
ate new links in the food web, perturbing the system dy- 
namics. When enough links are added, feedback loops will 
form, and the simulated ecology will suffer a mass extinc- 
tion. Over time, the system self organises to a state where 
the introduction of new species will be balanced by extinc- 
tions, and the system diversity fluctuates around some mean 
value. 

But what is this state that the system self organises to? 
The first suggestion was a critical state (Bak and Sneppen, 
1993), characterised by long range influences of a species 
extinction over others in the food web. The original model of 
Bak and Sneppen used to illustrate this idea is no more than a 
cartoon. The interactions between species in this model had 
no relation to biological interactions. The first attempt to 
use some real biologically inspired dynamics was probably 
Ecolab (Standish, 1994), which employed the well known 
Lotka-Volterra equations, for which a quite a bit of theoret- 
ical information is available. This model clearly self organ- 
ises to a state where speciation is balanced by extinction of 
average (Standish, 1999), although a variation of the model 
(incorporating a mechanism of specialisation) produces un- 
bounded growth in diversity (speciation exceeding extinc- 
tion) (Standish, 2002). 

So is this state a critical state? One problem is that criti- 
cality in self-organised systems is only achieved in the limit 



of zero driving rate — in this case zero mutation rate. Sole 
et al. (Sole et al., 2002) prefer the term self-organised in- 
stability. Whilst I am sympathetic to this notion, I would 
also like to point out that stability is very precise term in 
dynamical systems theory, referring to the behaviour of the 
linearised system around an equilibrium point. Unstable 
ecosystems do not have to fall apart — the classical Lotka- 
Volterra (Maynard Smith, 1974) limit cycle is a case in point. 
Rather the notion of an ecosystem persisting in time without 
falling apart is captured by permanence, for which a few 
modest results are known for Lotka-Volterra systems (Law 
and Blackford, 1992). So perhaps self-organised imperma- 
nence would be a more accurate description. 

Self organised critical systems are characterised by a 
power law distribution of extinction avalanches, and also 
a power law distribution of lifetimes. Traditionally, the 
presence of power law signatures in a self-organising sys- 
tem is taken as evidence of self-organised criticality. New- 
man (1997) developed another toy evolutionary model that 
exhibited power law spectra, with neither self-organisation 
nor criticality in sight. However, when the artificial con- 
stant diversity restriction is lifted in the obvious way, self- 
organisation reappears (Standish, 1999), and the model can 
also be understood as a mean field approximation of coevo- 
lutionary system that potentially admits critical behaviour. 

Ecolab demonstrates power law spectra of lifetimes (Stan- 
dish, 1999), with an exponent of -1. However, it has proven 
very difficult to measure the distribution of extinctions, as 
extinction avalanches overlap in Ecolab due to the finite rate 
of speciation. Conversely, studies of a similar model called 
Webworld claim an absence of any power law signatures 
(Drossel et al., 2001). I have implemented the Webworld 
model using the Ec<Lab (Standish, ) simulation system. I 
was similarly unable to see evidence of power law signa- 
tures, and propose a possible explanation. 

In this paper, I show that the Fourier transform of the di- 
versity time series is related to the lifetime distribution. Fur- 
thermore, in the limiting case of infinitesimal speciation, this 
transform is the distribution of extinction avalanches (ex- 
tinction frequency). 



Ecolab model 

We start with a generalised form of the Lotka-Volterra equa- 
tion 

ni = r i n i -n i '^2$ijnj. (1) 

j 

Here n, is the population of species i, r, is the difference 
between reproduction and death and p, 7 - is the interaction 
between species ;' and j. 

Periodically, each species i generates a number of mutant 
species, proportional to n,r,yu;, where /i, is the mutation rate 
for species i. For each mutant species, the parameters r„ 
Py , and fj( are mutated from the parent species according 
to additive or multiplicative processes — the exact details 
aren't important here, but are described in (Standish, 1994). 

One crucial property that is preserved by the mutation op- 
erator is boundedness (Standish, 2000). Boundedness en- 
sures that population sizes in eq (1) can never exceed a par- 
ticular limit. 

It turns out that a necessary condition for permanence in 
eq (1) is that the matrix p has positive determinant (Law and 
Blackford, 1992). The determinant can be written as a sum 

det|P|= Yl (-ir (/,) Pi Pl P 2; , 2 ---P n p„ (2) 

pGperm(l...,«) 

where perm(l ... ,n) is the set of permutations of the num- 
bers 1 . . . ,n and s(p) is the number of swaps involved in the 
permutation. 

All diagonal terms of p must be positive to ensure bound- 
edness of eq (1). 

Now consider permutations with one swap (z — » j,j — > ;'. 
If the terms p, 7 and p 7 , are of opposite sign (predator-prey 
case), then the contribution to the determinant is positive. 
However, if the terms have the same sign, (eg +ve, the mu- 
tual competition case, an increase in n,- causes « ; - to decrease, 
which reduces competition on n,, reinforcing the original 
change) then it describes a positive feedback loop between 
species i and j. 

Likewise, it can be seen that the term T = 
(— l)' s (p)p lpi p2 P2 • • • p„ Pn describes an s(p) feedback 
loop through the ecosystem, which is a negative feedback 
loop if T > 0, and a positive feedback loop if T < 0. The 
necessary condition for permanence can be interpreted as 
saying that negative feedback loops must dominate over 
positive feedback loops for the ecosystem to be permanent. 

As species are added to the system through speciation, 
new links are added to the foodweb at random. The chance 
of a positive feedback loop forming increases dramati- 
cally as the foodweb approaches its percolation threshold 
(Green and Klomp, 1999). Once this happens, an extinc- 
tion avalanche is almost certain. The twin pressures of spe- 
ciation and extinction through impermanence oppose each 
other leading to a state where the food web lies on its perco- 
lation threshold. Newth et al. (2002) examined the scaling 
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Figure 1: Lifetime distributions for different values of the 
maximum mutation rate /j in Ecolab. Histograms have been 
scaled so that the peaks are comparable in size. 

structure of the Ecolab model, and observed the critical be- 
haviour here. This is the strongest evidence yet that Ecolab 
self-organises to a critical state. 

Plots of the lifetime distribution for several different val- 
ues of the maximum mutation rate (mutation rates in Eco- 
lab are allowed to vary, but can never exceed the maximum 
value) are shown in figure 1 . These can be compared with 
other published data, such as (Standish, 1999). At higher 
mutation rates, the distributions exhibit a power law tail 
with exponent — 1 . As the mutation rate is turned down, the 
power law tail disappears, leaving a lognormal distribution. 
It is unclear whether the power law has disappeared alto- 
gether, or whether with the collection of more data it will be 
resolved out of the noise at the base of the graph. 

Relationship between diversity time series and 
lifetime distribution 

When a species becomes extinct, it may trigger secondary 
extinctions in other species, in a chain of extinctions known 
as an extinction avalanche. In the Bak-Sneppen model, these 
avalanches follow a power law distribution in avalanche size 
with exponent — 1 . However, it only becomes meaningful 
to discuss avalanche size in the limit of infinitesimal mu- 
tation rate, as otherwise the extinction avalanches overlap 
each other. In the Ecolab case, speciation occurs continu- 
ously, as do the resulting extinctions. More interesting is 
to discuss the frequency spectrum of extinctions, obtained 
by Fourier transforming the extinction time series. Diversity 
(number of species in the ecosystem at any point in time) 
is simply the difference between the speciation and extinc- 
tion time series — in the infinitesimal speciation limit, the 
diversity spectrum is identical to the extinction spectrum. 

The diversity time series can be written as a sum over spe- 
ciation events Sj and associated lifetimes if. 

D(t) = Y®(t-Sj)-®(t- Sj -Xj), (3) 

j 

where 
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Figure 2: Fourier transform of a typical diversity time series 
in Ecolab showing hyperbolic behaviour. 
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Figure 3: Fourier transform of a typical diversity time series 
in Webworld showing hyperbolic behaviour. 
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is the usual Heaviside step function. Taking the Fourier 
transform of this (and ignoring constant factors): 



D(k) = y^exp(z'fay) 



1 — exp(ikij) 
ik 



(4) 



Now we assume that speciation events occur every timestep 
(ie Sj = j), and that the lifetimes x ; are sampled from a nor- 
malised lifetime distribution T(z). Integrating over this dis- 
tribution yields: 



D(k) = J2^P(ikj) 



1 - f(k) 1 - T(k) 



ik 



(l-e ik )ik' 



(5) 



Now, if lirn^o T(x) < °°, then T(k) 
can predict that asymptotically, 

\D(k)\ - AT 1 ask- 



as k — > °°. So we 



(6) 



As k -> 0, 1 - e ik ik; 1 - f (Jfc) ~ (z)ik, where (x) is the 

mean of T(x). Even though a power law lifetime distribution 
would lead to an infinite (x), and any experiment, there is an 
upper cutoff to the lifetimes observed, which would reflect 
a finite (x). Therefore also D(k) ~ k^ 1 as k — > 0. Figures 
2 and 3 show D(k) for typical Ecolab and Webworld runs 
respectively, and both data sets demonstrate this hyperbolic 
law. The power law observed in the time series spectra is 
uninteresting, as it is a general feature of all such stochastic 
processes. Adami (1998) makes a similar point in his book 
— that power laws in the time series spectra are necessary, 
but not sufficient for self-organised criticality. 

Webworld 

The Webworld model was introduced by Caldarelli et al. 
(1998), with some modifications described in Drossel al. 
(2001). The model implemented here is taken verbatim from 
the latter paper, so I will give only a brief synopsis of the 



model here. The source code is available as part of the 
Ecq_ ab i 

software suite. 
Webworld has a population dynamics which is a generali- 
sation of the Lotka Volterra dynamics (eq 1) used in Ecolab: 



-n i + hi l 'Y J gij{t) -J2 n j8ji(t)- 



(7) 



This equation is called a functional response equation. X 
is a model parameter called ecological efficiency, and usu- 
ally taken to be X = 0.1. «o = R/X is a special species, 
called the environment. By choosing gij(t)ni — gji{t)rij = 
fiijninj, \/i,j > and g® = 1, equation (1) is recovered. 
However, unlike Ecolab, Webworld tracks resources, and so 
Hi is perhaps better interpreted as the amount of biomass rep- 
resented by species i than a population size. 

In Webworld, the functional response term gij is given by 



Sij{t) 



S u f u (t)nj(t) 



bnj{t) + J2 k UkiSkjfkj (t)n k (t) 
where the efforts fij are given recursively: 



fij(t) 



suit) 



(8) 



(9) 



Drossel et al. (2001) show that allowing species to vary the 
amount of effort in this way is an evolutionary stable strat- 
egy. The 0Cy < 1 terms above represent that different species 
do not compete as strongly as members of the same species 
(Ofl = l, Vi). 



-c+(l -c)q t 



(10) 



where < c < 1 is a competition parameter that strongly 
influences the final steady state diversity of the model. The 
precise definition of the interaction terms 5,/ and qij is very 
interesting, but not germane to the argument here. 

In (Drossel et al., 2001), the equations are evolved in time 
until the ecosystem reaches equilibrium, or until a large pe- 
riod of time has elapsed before another species is added to 

1 http://parallel.hpc. unsw.edu.au/ecolab 
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Figure 4: Lifetime distributions for different values of c for 

rr/e^syffiJfii. This is to mimic the "infinitesimal speciation 
rate" mentioned in relation the self-organised criticality. In 
this study, a single species is added periodically every 20 
time units, which was chosen empirically as being suffi- 
ciently rare for the ecosystem to stabilise between specia- 
tions. 

Figure 3 shows the spectrum of the diversity time series 
for a run with c = 0.4. Figure 4 shows a lifetime distribution 
from these runs, illustrating that the best functional form is a 
lognormal distribution of lifetimes. Its important to note that 
lognormal distributions are often confused with power law 
distributions (Mitzenmacher, 2003), however are distinctly 
different in that lognormal distributions have a mean (x), 
whereas power law distributions do not. However, as men- 
tioned in the section describing the Ecolab model, very low 
speciation rates were used in this model, and it is possible 
that a power law tail is hidden within the noise at tail of the 
graph. 

Quite another explanation for the absence of critical be- 
haviour in Webworld comes from considering the dynam- 
ics of the effort coefficients fy. These are iterated within a 
timestep to determine evolutionary stable values — this is a 
model of how predators exhibiting high phenotypic plastic- 
ity might adapt to variations in food sources. It would seem 
to be a less plausible model of less complex organisms that 
might not be so choosy about their food source. Whatever 
the biological realism of this process, the fij coefficients de- 
scribe the effective foodweb, which tends to be quite sparse 
and without many loops (Quince et al., 2002). Could it be 
that this process prevents the percolation threshold of the 
foodweb from being reached? 

Mean Field Ecolab model 

In (Standish, 1999), I introduced a mean field 2 version of 
the Ecolab model (which I dubbed Ecolab — ). That model 
is a simple multiplicative process, which is related by log- 



2 Some people use mean field to mean panmictic. The Ecolab 
model described in previously is already panmictic. By mean field, 
I mean that each species experiences a stochastic force that is the 
average of the interspecific interactions in the full model. 



arithms to the standard isotropic ID random walk process. 
The lifetime distribution is known as the first-passage time 
distribution in this subject, and is known to exhibit a x~ 3 / 2 
tail (Redner, 2001). Note that this is different from, but still 
compatible with, the upper bound of x~ 1 derived in (Stan- 
dish, 1999). 

However, this model has a lognormal limiting distribution 
of population sizes n of the form: 



. . 1 . (lnn-rf) 2 . 

P{n,t) = —F^M t: ) 

nyjt At 



(11) 



This distribution does not satisfy boundedness. 

In order to introduce boundedness, we need to reintroduce 
the quadratic term into the mean field model: 



h = n(r — pn +y), 



(12) 



where y is an uncorrected stochastic variable, with zero 
mean. 

Taking logarithms t, = lnn and applying the Ito transfor- 
mation formula (Karlin and Taylor, 1981, p. 372, e.g.), eq 
(12) can be written as a Langevin equation: 



K = (r-\-^ + y). 



(13) 



The extra term of \ comes from the effect of change of 
variables on the stochastic term y. Langevin equations can 
be converted into a Fokker-Planck equation describing the 
probability distribution w(x,t) that ^ has the value x at time 
t (Risken, 1984): 



dw d 2 w d , 1 _ r . 
^ = W-Jt {r -2-^ e )W 



(14) 



Taking the Laplace transform of equation (14) yields a 
second order homogeneous ordinary linear differential equa- 
tion: 

|- (w'(x,s) - (r- l - - pe>) +sw = 0. (15) 

The full time dependent equation doesn't appear to be 
amenable to analytic treatment, however the time indepen- 
dent equation (s = 0) can be reduced to a 1st order ODE. 
Let 



v(.v) = exp( (r-^)x-$e } 



(16) 



r---W)y{x)=g(x)y(x), (17) 



and write wq(x) = w(x,0) = y(x)v(x). Substitute this into 
equation (15), and one obtains: 

d 



dx 



iyv') 



o 



, f dx ' * 



A {-$) r -±r{ l --r-$e x )+A x (18) 



where F(a,x) is an incomplete gamma function 
(Abramowitz and Stegun, 1965, 6.5.3), and Aj are constants 
of integration. Substituting this into the expression for 
p(n,0) yields: 

p(n,0) = -w(lnrc,0) = 
n 

„H e -P»^Ai+Ao(-pr'r(l-r,-P«)) (19) 



From the series F(a,x) ~ F(a) 
Ryzhik, 1980, 8.354), one can see: 



■ ■ ■ (Gradsteyn and 



p(n,0) ~ n r -h-V"(A (-$)'- ?r( 1 --r)+A^j + 



A 



n(r-i) 



(20) 



which is normalisable if and only if r > j and Aq = 0. We 
may therefore set wq(x) = y(x). 

The asymptotic behaviour at large times translates into the 



the small s regime. We can compute w\ (x) — — 
differentiating eq (15) with respect to s. 



s=0 



by 



A(w' 1 -(r-i-P^)w 1 )+w = (21) 



to which the solution is: 



w 1 (x)=w (x)(a 1 +J ^-(x')J w (x")dx>^J. (22) 

The innermost integral can be evaluated, the answer being 
(Abramowitz and Stegun, 1965, 6.5.2) 



w (x')dx' = r r+ h(r-l,^) (23) 



with y being another of the incomplete gamma functions. 
This can be represented as a series (Gradsteyn and Ryzhik, 
1980): 



/ 



^»>^ r -ii: '-')yr- (24) 

^ n\{r-l+n) 
Performing the integral on each term of the series yields: 

wi(x)=w (x)x (25) 

where E\(x) is the exponential integral (Abramowitz and 
Stegun, 1965,5.1.1). 

Interestingly, for the special case r — j e Z + , the result 
can be expressed as a finite series. It might seem that r can 



be chosen to any value by scaling the time dimension with- 
out loss of generality, however that is not the case, as the 
timescale is already set by the variance of the stochastic term 
yin eq. (12). 

Considering the special case r = |, and making use of the 
identity y(l,x) = (1 - e~ x ) (Gradsteyn and Ryzhik, 1980, 
8.352), we have: 

wi (x) = wo (x) (A i - p- l e- x - r( - 1 , pe x ) ) (26) 

To compute the asymptotic form of the lifetime distribu- 
tion, we make use of first passage theory (Redner, 2001). 
The first passage probability F(n,t\no) of the population 
having value n, given a starting value no at t = is related to 

p(n,t\n ): 

p(n,t) = 8(f )8(n — no) + / F(n,t — t'\no)p(no,t\no)dt' 
Jo 

(27) 

Taking the Laplace transform, and rearranging gives us 



F(n,s) = p(n,s)/p(n ,s) n Q ^n 



(28) 



(as we're not interested in the n = no case). For concreteness, 
let no = 10 and n = 1, as is taken in the case of Ecolab ex- 
periments computing the lifetime distribution (Fig. 1). The 
asymptotic form can be computed directly from F(n,s): 



? («,1/T) = J 



e~ zt F(n,t)dt ■ 



• J F(n,t)dt 



F(«,t)~-F(«,1/t) 
p(no,0)^£l 



(29) 



(30) 



s=0 



-p(n,0)^ 



s=0 



1 2 p(n Q ,0) 2 

Unless the numerator vanishes, the long time tail will obey 
a x~ 2 power law. Substituting equations (16) and (26) yields 
for the case r= |: 



F(n,x) 



,2p« 



^-^-+r(-i,» )-r(-i,») 

P«o pn 



(31) 

Since T(— l,n) > r(— l,«o) f° r n < n o, this derivative term 
is negative. It seems unlikely to vanish for any value of r > 
0, however this will need to be checked numerically. 

This result is interesting. The mean field model can be 
considered as a neutral model, in the sense of the neutral 
shadow models proposed by Bedau and Packard (Bedau 
et al., 1998). An observed excess of lifetimes over the mean 
field case (in Fig lar 1 distribution is observed) would rep- 
resent coadaption by the species in the ecosystem. 

Conclusion 

In this paper, I consider the question of self-organised criti- 
cality in a couple of evolutionary ecology models (Ecolab 



and Webworld). In spite of their similarity, only Ecolab 
appears to self-organise to criticality, whereas Webworld's 
self-organised state appears to be noncritical, in agreement 
with Webworld's creator's statements. 

Whilst it is possible that experiments have not been run 
long enough to observe critical behaviour, a more likely ex- 
planation is organismal plasticity in Webworld prevents long 
range interdependence of species in the foodweb from build- 
ing up. 

A mean free approximation to the Ecolab model is solved, 
and the lifetime distribution from this model is expected to 
have a T~ 2 asymptotic behaviour. The fact that the real Eco- 
lab model appears to have a T~ 1 asymptotic behaviour hints 
at adaption occurring within that system. 

Finally, the spectral density of the diversity time series 
(which is related to the distribution of extinction avalanches) 
is expected to have a 1/f behaviour, regardless of the un- 
derlying process, so this should not be taken as evidence of 
self-organised criticality. 
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